      !< @file
      !> @brief Subrotina para o cálculo do erro nas medidas referente ao modelo aproximado obtido com a inversão.
      !>
      !> Calcula o erro dos sismogramas referentes ao modelo obtido com a inversão em relação aos sismogramas dados como entrada para a inversão.

      subroutine calcMisfitFinal(gruposSR,Nsource,Nrec,Nttotal,ntps,unit,freqmax)

      implicit none


      double precision, parameter :: Pi = 4.*atan(1.) !< número Pi


      integer :: gruposSR !< Número de grupos de pares fonte/receptores
      integer :: Nsource(gruposSR) !< Número de fontes em cada grupo de pares fonte/receptor
      integer :: Nrec(gruposSR) !< Número de receptores em cada grupo de pares fonte/receptor
      integer :: Nttotal !< Número de instantes de tempo gravados
      integer :: ntps !< amostragem, número de instantes de tempo por segundo
      integer :: unit !< unidade onde serão impressos os resíduos
      double precision :: freqmax !< Frequência máxima em Hz para a qual será calculado o erro

      real, dimension(:,:), allocatable :: Pt !< sismogramas (traços) de cada arquivo, referente ao modelo aproximado encontrado na inversão
      real, dimension(:,:), allocatable :: Ptbarra !< sismogramas (traços) de cada arquivo, referente aos dados da inversão  
      double precision, dimension(:,:), allocatable :: dPt !< conversão de Pt em dupla precisão
      double precision, dimension(:,:), allocatable :: dPtbarra!< conversão de Ptbarra em dupla precisão


      integer :: igsr !< contador de grupos de pares fonte/receptor
      integer :: ierr !< indicador de erro na alocação
      integer :: ks !< contador de fontes
      integer :: kr !< contador de receptores
      integer :: it !< contador de intantes de tempo
      integer :: kw !< contador de frequências

      double precision:: MisfitT !< Resíduo entre os sismogramas medidos e calculados (dados no tempo) 
      double precision:: aux
      double precision, dimension(:), allocatable:: MisfitW !< Resíduo entre os sismogramas medidos e calculados (dados na frequência) 
      double precision, dimension(:), allocatable:: auxW 

      complex*16, dimension(:), allocatable:: Pwaux !<transformada de um traço dos sismogramas do modelo aproximado para o dominio da frequencia
      complex*16, dimension(:), allocatable:: Pwauxbarra !<transformada de um traço dos sismogramas do modelo exato para o dominio da frequencia

      integer :: nomega !< número de frequências para as quais serão calculados os resíduos
      double precision, dimension(:), allocatable :: Omegas !< frequências angulares


      nomega = int(freqmax)

      allocate(omegas(Nomega))!frequências angulares para as quais serão calculados os erros
      do kw=1,Nomega
	  omegas(kw) = dfloat(kw)*2.0d0*Pi
      enddo
      allocate(Pwauxbarra(Nomega))
      allocate(     Pwaux(Nomega))
      allocate(MisfitW(Nomega))
      allocate(auxW(Nomega))

      do kw=1,Nomega
	  MisfitW(kw) = 0.0
	  auxW(kw) = 0.0
      enddo



      ! Descompactar e ler sismogramas
      call system('tar -xzf sismD.tgz') 
      write(*,*)'sismogramas  da inversao descompactados'

      call system('tar -xzf sismA.tgz') 
      write(*,*)'sismogramas  do modelo aproximado descompactados'



      write(*,*)'Calculando erros...'

      MisfitT=0.0
      aux = 0.0
      do igsr=1,gruposSR ! Loop no número de grupos fonte/receptor

	  allocate(Ptbarra(nttotal,nrec(igsr)),stat=ierr)
	  if(ierr.ne.0)then
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif

	  allocate(Pt(nttotal,nrec(igsr)),stat=ierr)
	  if(ierr.ne.0)then
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif

	  allocate(dPtbarra(nttotal,nrec(igsr)),stat=ierr)
	  if(ierr.ne.0)then
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif

	  allocate(dPt(nttotal,nrec(igsr)),stat=ierr)
	  if(ierr.ne.0)then
	      write(*,*)'.........allocation error',ierr     
	      stop
	  endif

	  do ks=1,Nsource(igsr) ! Loop nas fontes
	      call loadSismog(Ptbarra,nttotal,nrec(igsr),igsr,ks,'sismD')
	      call loadSismog(Pt     ,nttotal,nrec(igsr),igsr,ks,'sismA')

! 	      open(15,file='lixo.pld')
! 	      do it=1,nttotal
! 		  write(15,*)Ptbarra(it,15),Pt(it,15)
! 	      enddo
! 
! 	      close(15)
! 	      call system('gnuplot lixo.plt')
! 	      read(*,*)


	      do kr=1,Nrec(igsr)
		  do it=1,nttotal
		      MisfitT = MisfitT + ((Ptbarra(it,kr)-Pt(it,kr)))**2
		      aux = aux + (Ptbarra(it,kr))**2
		  enddo
	      enddo


	      do kr=1,nrec(igsr)
		  do it=1,nttotal
		      dPtbarra(it,kr) = DBLE(Ptbarra(it,kr))! isso pq a rotina dft foi feita para double precision
		      dPt(it,kr)      = DBLE(Pt(it,kr))! isso pq a rotina dft foi feita para double precision
		  enddo
	      enddo


	      do kr=1,Nrec(igsr)
		  call dft(Pwauxbarra,omegas,Nomega,ntps,dPtbarra(1,kr),nttotal)
		  call dft(Pwaux     ,omegas,Nomega,ntps,     dPt(1,kr),nttotal)
		  do kw=1,Nomega
		      MisfitW(kw) = MisfitW(kw) + DCONJG(Pwauxbarra(kw)-Pwaux(kw))*(Pwauxbarra(kw)-Pwaux(kw))
		      auxW(kw) = auxW(kw) + DCONJG(Pwauxbarra(kw))*(Pwauxbarra(kw))
		  enddo
	      enddo


	  enddo! fim do loop nas fontes

	  deallocate(Ptbarra)
	  deallocate(Pt)

      enddo ! fim do loop no número de grupos fonte/receptor

      MisfitT = MisfitT / aux
      write(unit,*) '#Resíduo com os dados no domínio do tempo: ',MisfitT,sqrt(MisfitT)*100.0,'%'

      write(unit,*) '#Resíduo na frequencia--  freq em Hz versus resíduo, erro relativo em %'
      do kw = 1,Nomega
	  MisfitW(kw) = MisfitW(kw)/auxW(kw) 
	  write(unit,*) omegas(kw)/(2.0d0*Pi),  MisfitW(kw), sqrt(MisfitW(kw))*100.0     
      enddo




      deallocate(omegas)
      deallocate(Pwauxbarra)
      deallocate(Pwaux)
      deallocate(MisfitW)
      deallocate(auxW)

      return
      end subroutine calcMisfitFinal